A Two-dimensional Superconductor in a Tilted Magnetic Field - new states with 

finite Cooper-pair momentum 
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Varying the angle between applied field and the conducting planes of a layered superconductor 
in a small interval close to the plane-parallel field direction, a large number of superconducting states 
with unusual properties may be produced. For these states, the pair breaking effect of the magnetic 
field affects both the orbital and the spin degree of freedom. This leads to pair wave functions 
with finite momentum, which are labeled by Landau quantum numbers < n < 00. The stable 
order parameter structure and magnetic field distribution for these states is found by minimizing the 
quasiclassical free energy near Hc2 including nonlinear terms. One finds states with coexisting line- 
like and point-like order parameter zeros and states with coexisting vortices and antivortices. The 
magnetic response may be diamagnetic or paramagnetic depending on the position within the unit 
cell. The structure of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states at = is reconsidered. 
The transition n ^ 00 of the paramagnetic vortex states to the FFLO-limit is analyzed and the 
physical reason for the occupation of higher Landau levels is pointed out. 
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I. INTRODUCTION 

In this paper a theoretical study of a two-dimensional, 
clean-limit superconductor in a tilted magnetic field is 
presented. Such systems exist in nature; several classes of 
layered superconductors of high purity with conducting 
planes of atomic thickness and nearly perfect decoupling 
of adjacent planes have been investigated in recent years. 
These include, among many others, the intercalated tran- 
sition metal dichalcogenide TaS2 — {pyridine), the or- 
ganic superconductor n - {BEDT - TTF)2Cu{NCS)2, 
and the magnetic field induced superconductor A — 
{BETS)2FeCh. 

Depending on the angle 6 between applied field and 
conducting planes the nature of the pair-breaking mech- 
anism limiting the superconducting state can be contin- 
uously varied. For large 6 the usual orbital pair-breaking 
mechanism dominates and the equilibrium state is the 
ordinary vortex lattice. With decreasing 0, in a small 
interval close to the parallel direction, spin pair-breaking 
becomes of a magnitude comparable to the orbital effect 
and both mechanisms must be taken into account. For 
the plane-parallel field direction, = 0, the orbital ef- 
fect vanishes completely and the superconducting state 
is solely limited by paramagnetic pair breaking. The su- 
perconducting state expected in this limit is the Fulde- 
Ferrell-Larkin-Ovchinnikov (FFLO) state^-^ The tilted- 
field arrangement, which allows to control externally the 
relative strength of both pair-breaking mechanisms, has 
first been investigated by Bulaevskii'^. 

The upper critical field Hc2, where a second order 
phase transition between the normal-conducting and 
the superconducting state takes place, has been calcu- 
lated for arbitrary angle 9 and temperature T = by 
Bulaevskii'^. This treatment was generalized to arbitrary 
T by Shimahara and Rainer^ . The field Hc2 has a cusp- 



like shape, considered both as a function of or T, with 
different pieces of the curve belonging to different val- 
ues of the Landau quantum number n (n = 0, 1, . . .). In 
the orbital pair breaking regime, for large 6, one finds 
as expected n = 0. As is well known, this lowest value 
n — Q determines the (orbital) upper critical field of the 
familiar vortex state, both in the framework of Ginzburg- 
Landau(GL)- and microscopic theories of superconduc- 
tivity. With decreasing 9, higher-n segments of the criti- 
cal field curve appear close to the plane-parallel orienta- 
tion. For — + one finds'* n ^ 00 and agreement with 
the FFLO upper critical field. Thus, in this purely para- 
magnetic limit, the stable state below Hc2 must be the 
FFLO state. 

Paramagnetically-limited superconductivity differs in 
fundamental aspects, such as Meissner effect and spin- 
polarization, from the behavior of the usual, orbitally- 
limited superconducting state. In the FFLO state pair- 
ing takes place between electrons with momentum and 
spin values (fc + (f/2,|) and (— fc + q/2,l). This leads 
to Cooper pairs with finite momentum Tiq and a spa- 
tially inhomogeneous superconducting order parameter 
given by A(r) = Aq exp(«gr^ (or by linear combinations 
of such terms with the same absolute value of q). The 
pair-breaking is entirely due to the Zeeman coupling be- 
tween the magnetic moment jj. of the electrons and the 
external magnetic field H. The general rule, for bulk 
superconducting states, that gradient terms in the free 
energy must only be taken into account if a nontrivial 
vector potential is present breaks down for the FFLO 
state. 

At T = 0, the Cooper pair momentum of the FFLO 
state is approximately given by fiq = \pf'\ ~ PFi\, where 
PFi \ ^ fJ-H y^2m/ Ep is the difference in Fermi mo- 
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mentum between spin-up and spin-down electrons. With 
increasing T the FFLO wave number q decreases and van- 
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ishcs at the tricritical point Tt„ = 0.56 T^. The FFLO 
state is only stable for T < T^j-i, where its upper criti- 
cal field Hp FLO exceeds the Pauli limiting field Hp of 
the homogeneous superconducting state^'^. At T = 0, 
IJ,Hp = Aq/V^, where Aq is the superconducting gap at 
T = 0. The second order phase transition line HppLoiT) 
depends on the shape of the Fermi surface. In this paper 
we use a cylindrical Fermi surface appropriate for a two- 
dimensional (2D) geometry. The corresponding critical 
field^ is given by ^HppLo = Aq at T = 0. 

Between the ordinary vortex state with n = and the 
FFLO state with n ^ oo a countable infinite number 
of unconventional superconducting states, characterized 
by Landau quantum numbers n = 1,2,..., exist. The 
transition from the vortex state to the first of these, the 
n = 1 state, occurs at an angle 9i given approximately 
by 

sm-Oi^-r—K (1) 

Hp mvp 

where H°2^ and Hp arc the "pure" orbital and param- 
agnetic upper critical fields, respectively. Since Hp ^ 
H°2^ the experimental upper critical field for a three- 
dimensional sample is given by H°2^. Because 6*1 ^ 1 
(generally 0i will be of the order of magnitude of 1 De- 
gree), the perpendicular component H± = HsinO for all 
of these states with n > will be much smaller than 
the parallel component = H cos 9. Thus, these states 
will have some properties in common with the FFLO 
state, namely strong paramagnetic-pair breaking, a spa- 
tially inhomogeneous order parameter, and Cooper pairs 
with finite velocity of the center of mass coordinate. De- 
spite this similarity with regard to general features, the 
order parameter structure for the n > states may be 
completely different, even for large n, from the FFLO 
state. The reason is, that a finite perpendicular compo- 
nent H±, no matter how small, implies a new and rather 
stringent topological constraint on the equilibrium struc- 
ture, namely the flux quantization condition. The sub- 
ject of the present paper is the detailed investigation of 
the structure of these n > states, which might be re- 
ferred to either as FFLO precursor states or as param- 
agnetic vortex states, in the vicinity of the upper critical 
field Hc2- A theoretical treatment of these FFLO precur- 
sor states, reporting several essential results and an out- 
line of the calculation, has been published previously*. 
This paper* will be referred to as KRS in what follows. 
In the present paper many new results are reported and 
the treatment is extended with regard to several points, 
including finite values of k, the purely paramagnetic limit 
= 0, and the transition n — !■ oo. 

It should be pointed out, that the physical origin of the 
Landau level quantization effects for Cooper pairs, con- 
sidered in the present paper, is very different from the 
Landau quantization effects for single electron states dis- 
cussed in a large number of publications by Tesanovic et 
al.^, Rajagopal et al.^", Norman et al.'^-'^ and others. The 
latter are mainly concerned with the relative-coordinate 



degree of freedom of the two bound electrons constitut- 
ing a Cooper pair and lead to measurable consequences 
only outside the range of validity of the quasiclassical ap- 
proximation, at very low temperature T < {kpTc)^ /Ep 
and/or high fields. In addition, a mechanism is re- 
quired to suppress the Zeeman effect, which is neglected 
in the theoretical treatment and is not compatible with 
the predicted phenomena. The question whether the 
most dramatic consequences^^ (reentrant superconduc- 
tivity) of this type of Landau quantization effects will 
be observable, has been the subject of a controversial 
discussion^'"*'^^. In contrast, the present Landau level 
quantization mechanism is a consequence of the Zeeman 
effect, concerns the center of mass motion of the Cooper 
pairs, and can be described (as will be discussed shortly) 
by means of the quasiclassical theory of superconductiv- 
ity. 

Restricting ourselves to the vicinity of the upper criti- 
cal field Hc2 we may use an expansion of the free energy 
in powers of the order parameter A, keeping only a fi- 
nite number of terms. An analogous gradient-expansion, 
which would lead to a relatively simple GL-like theory 
with a finite number of spatial derivatives of A, does, 
unfortunately, not exist for the present problem. Such 
an expansion may be performed for ^ = 0, in the purely 
paramagnetic limit, near the tricritical point Tiri, where 
the order parameter gradient is small because the char- 
acteristic length of the FFLO state diverges at Ttri- 
However, for finite H± a small characteristic length for 
order parameter variations does not exist in the relevant 
range of temperatures, and the spatial variation of A 
must be taken into account exactly. One might still 
hope that a GL theory with a finite number of deriva- 
tives, although not accurate, will be useful to predict the 
qualitative behavior of the superconducting states near 
Hc2 correctly; bearing in mind for example the results of 
standard GL for type II superconductivity. However, for 
the mixed orbital-paramagnetic pair-breaking phenom- 
ena under discussion, there is not even a single point 
on the temperature scale where a GL theory with a finite 
number of derivatives is valid. Such a theory is only valid 
near Tc where no FFLO state exists, or near Ttri in the 
"vicinity" of the paramagnetic limit, i.e. for extremely 
large n. The latter region is inaccessible both from a 
numerical and a experimental point of view. In this con- 
text, it should also be noted that the final equilibrium 
structures do not show any continuity with regard to n. 

Fortunately, the present problem does not require solv- 
ing the full set of Gorkov's equations because the simpler 
set of quasiclassical equations may be used instead, as 
pointed out by Bulaevskii"^. The large parallel compo- 
nent of the applied magnetic field, acting only on the 
spins of the electrons, is exactly taken into account by the 
Zeeman term. Thus, with regard to this component no 
question, as to the validity of the quasiclassical approxi- 
mation, arises. The magnitude of the perpendicular com- 
ponent H±, on the other hand, must obey the usual qua- 
siclassical condition hujc < kpT, where ojc = eH±/mc, 
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or sm6{ehH / mc) < ksT. Inserting the highest possible 
field H = Hp in the latter relation, one finds that the 
quasiclassical approximation holds indeed for not too low 
temperatures, T/Tc > ksTc/Ep, in the interesting range 
of tilt-angles < 0\, where the new paramagnetic vortex 
states appear. 

In most papers on paramagnetic pair-breaking and the 
FFLO state the influence of orbital pair-breaking is com- 
pletely neglected. This means, that the GL-parameter k 
tends to infinity and that all spatial variations of the mag- 
netic field can be neglected. For three-dimensional super- 
conductors this approximation implies that the orbital 
critical field is much higher than the paramagnetic Pauli- 
limiting field. This is impossible to achieve^^ for BCS-like 
superconductors, because the superconducing coherence 
length cannot be smaller than an atomic distance. It 
seems unlikely even for unconventional materials^® where 
many-body effects may lead to a strong renormalization 
of the input parameters. For the present 2D situation, 
the suppression of the orbital pair-breaking effect is en- 
tirely due to geometrical reasons, and no restriction on 
the value of k is required in order to reach the purely 
paramagnetic limit at parallel fields. Thus, keeping all 
terms in the quasiclassical free energy related to spatial 
variations of the magnetic field, will allow us to study 
type II superconductors with arbitrary n or even type I 
material. Large-K superconductors show, however, still a 
practical advantage because of their larger critical angle 
6*1 [see Eq. (1)]. 

This paper is organized as follows. In section II Eilen- 
bergers quasiclassical equations generalized with regard 
to a Zeeman coupling term, as well as the corresponding 
free energy functional, are reported. The expansion of 
the free energy near the upper critical field, for a general 
2D quasi-periodic state, is treated in section III. Two 
limiting cases of the analytical results, the GL limit and 
the structure of the ordinary vortex lattice, are reported 
in appendices. The numerical results for the paramag- 
netic vortex states, at finite perpendicular field, are re- 
ported and discussed in section IV. The structure of the 
FFLO state, for the special case of vanishing perpen- 
dicular field, is reconsidered in the present quasiclassical 
framework in section V. The non-trivial transition ^ ^ 
(or n ^ c») to the purely paramagnetic limit is analyzed 
in section VI. An explanation for the increase in n, in 
terms of the finite momentum of the Cooper pairs in the 
paramagnetic vortex states, is also reported in this sec- 
tion. The results are summarized in the final section VII. 



II. QUASICLASSICAL EQUATIONS WITH 
ZEEMAN TERM 



We need a weak- coupling, clean-limit version of the 
quasiclassical theory^^'^^, which contains all terms re- 
lated to the coupling of the electron's spins to an exter- 
nal magnetic field. A general quasiclassical theory which 
covers Zeeman coupling has been published by Alexander 



et.al'^^. The 4x4 Green' s function matrix appearing in 
this work may be considerably simplified for the present 
situation. Since we neglect spin-orbit coupling, the direc- 
tion of the magnetic induction B in spin-space may be 
chosen independently from the direction of B in ordinary 
space; we adopt the usual choice of B being parallel to 
the z— direction in spin space. Then, only six essential 
Green's functions remain, which are denoted by 



f + h 



/(+) 



9i+) 



9 + 93 



'(+) 
9{-) 



r + ft 

9-93- 



Here, /, /+, g denote the Greens functions in the absence 
of Zeeman coupling, and /a, ft, gs are the additional 
Green's function components in the the ^—direction of 
spin space. The three equations for the left group 
/(_), /(!^), g(-) are decoupled from the three equations 
for the right group /(+), 9(+) ^-nd differ only by 

a negative sign in front of the magnetic moment ji = 
h\e\/{2mc) of the electron. Also, for each group a sep- 
arate normalization condition g^^j + /(+)/('^-) = 1 and 

.J + /(_)/^^-) = 1 respectively, exists. Therefore, it 

is convenient to introduce Greensfunctions /, /+, g, de- 
fined by 

f+{f,k,u!s) = /(+)(f, fc,a;), 

g{f,k,iJs) = 5(_)(f, fc,a;), 

which are functions of the spatial variable f, the quasi- 
particle wave-number k, and the complex variable = 
CO + ifjiB. The 2D variable r denotes positions in the con- 
ducting (x, y)-plane. The real variable u) takes the values 
of the Matsubara frequencies oji = {2l+l)TTkBT; the Mat- 
subara index I will not always be written down explicitly. 
The second group of Green's functions /(+), f^L): 9{+) 
may be expressed by similar relations in terms of f,f~^,g 
if Us is replaced by lo*. 

Using the Green's functions /, /+, g, the quasiclassical 
equations with Zeeman coupling become formally similar 
to the quasiclassical equations without spin terms. The 
nonlinear transport equations for /, /+ are given by 



2ws + %VF{k)'dr f{f, k, (Js) = 2A{r)g{f, k, w^), 
2us - MfWCI f^ir,k,tOs) = 2A*{f}gif,k,uJs), 



(2) 



where the Green's function g is given by the normal- 
ization condition 



g{f,k,ujs) = 1 - f{r,k,ujs)f^{r,k,ujs) 



1/2 



(3) 



Here, vp{k) denotes the Fermi velocity and -dr is 
the gauge-invariant derivative defined by -dr = — 
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i{2e/hc)A. The order parameter A and the vector po- 
tential A must be determined selfconsistently. 



The self-consistency equation for A is given by 



J 



Nd 



2wkBT - + In (T/Te) J A(r) = tt/cbT Y^jSk' [/(f, fc', w,) + /(f, fc', a;,*) 



(4) 



r 



where Nd is the cutoff index for the Matsubara sums. 
The self-consistency equation for A is Maxwell's equation 



Vr X (-B(r) + 47rM(f)) = 



c ^ 

z=o 



4-77 



VFik')^g{f,k',uJs), 



(5) 



where TV^? is the normal-state density of states at the 
Fermi level. The r.h.s. of Eq. (5) is the familiar (orbital) 
London screening current while the magnetization M is 
a consequence of the magnetic moments of the electrons 
and is given by 



M(r) =2n'^NFB{r) 



Nd 



(6) 



The first term on the r.h.s. of Eq. (6) is the normal 



state spin polarization. The second term is is a spin 
polarization due to quasiparticles in the superconducting 
state. 

The following symmetry relations hold for solutions of 
Eqs. (2,3,4,5) 



g*if,-k,uj*) =g{f,k,ujs), 

f~^{r,k,uJs) = f*{f,~k,uj*), 

g{r, -k, -ujs) = -g{r, k, ujs), (7) 

/(f, -fc, -Us) = fir, k, ujs), 

f+{r,-k,-u;s) = /+(r,fc,a;s). 



which have been extensively used in the calculations de- 
scribed in the next sections. 

The quasiclassical equations (2,4,5) may be derived as 
Euler-Lagrange equations of the Gibbs free energy func- 
tional G, which is given by 



G 



TrksT U| +^^(^/^«)) \^\^-^kBTNF ^ j ^I{r,k,u;,) 

l=—oo / l= — oo 



r 



(8) 



The area of the sample is denoted by Fp and the 
fc— dependent quantity / is given by 

I{r, k, cOs) = A/+ + A*/ +(g-^^. 





n.vF^\ 







Mf 



An important reference state for the present problem 
is the purely paramagnetically limited homogeneous su- 
perconducting state, which is realized for our 2D super- 
conductor if the magnetic field is exactly parallel to the 
conducting planes. In this case, the vector potential and 



the gradient terms in the transport equations may be 
omitted. At T = the free energy difference between 
the superconducting and normal-conducting states may 
be derived analytically. It is given by 



Gs - Gn = Nf iti^H^ - Al/2) 



(9) 



and vanishes at the Pauli critical field Hp. For higher 
T the self-consistency equation for the gap must be 
solved numerically, yielding agreement with previous 
results. Let us investigate the magnetic response in 
this purely paramagnetic limit. It is neglected in most 
theoretical treatments, but is of particular interest if the 
influence of finite values of the GL-parameter k is to be 
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taken into account. To obtain the magnetization due 
to the spins, the coupled self-consistency equations (4,5) 
have to be solved. Using dimensionless quantities defined 
in appendix A the gap equation takes the form 



1 



while Maxwell's equation reduces to 




= 0, 
(10) 



Nd 



u)i + IIJ.B 



=0 + {uji+niBf 



(11) 

Note that the orbital screening current [the r.h.s. of 
Eq. (5)] is completely absent for the plane-parallel field 
direction. At T = the r.h.s. of Eq. (11) vanishes ex- 
actly. This means that the normal state spin polarization 
[first term on the r.h.s. of Eq. (11] is exactly canceled by 
the spin polarization due to the superconducting quasi- 
particles [second term on the r.h.s. of Eq. (11)]. The 
numerical solution shows that the quasiparticlc polariza- 
tion decreases with increasing T and vanishes at A = 0, 
where the magnetic behavior of the normal- conducting 
state is recovered. 

In the rest of this paper dimensionless quantities as 
introduced by Eilcnbcrger will be used. These quantities 
are listed in appendix A. Any exceptions will be men- 
tioned explicitly. 

In the next sections the stable order parameter struc- 
ture of a 2D superconductor in the vicinity of the phase 
boundary will be investigated. The phase boundary 
Hc2{T) itself is given by the highest solution of the 
equation'' 



= \ut+t 



f 

Jo 



ds 



1-e- 



1 



sinh st 

- cos{nHs)e-"^''/^L„{H±s'^/2) 



(12) 



where the integer n = 0,1,2,... is Landau's quantum 
number, lod is the Debye frequency, and L„ is a Laguerre 
polynomiaP^ of order n. A typical phase boundary is 
shown in Fig. 1. Each piece of the nonmonotonic Hc2 
curve is characterized by a single value of n. An infinite 
number of eigenstates 4>n,k exists, belonging all to the 
same, highly degenerate eigenvalue n. For the present 
gauge, these are given by 



(-1)" . _!l±(y 



(t>n,k{r) =A ^ e'—e 



k 
Hi 



(13) 



where fc is a real number and Hen is a Hermite 
polynomials^ of order n. The functions (13) are orthog- 
onal and normalized, 



{<l>n,k->4>m,l) = 5n,m5{k - I), 



(14) 
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FIG. 1: Phase boundary of the superconducting state at 
t = 0.1 for tilt angles © between 0.1 and 2.0 using a value 
jj, = 0.04 for the dimensionless magnetic moment of the elec- 
tron. The numbers 0, 1, 2, . . . are Landau quantum numbers 
characterizing the individual pieces of the curve. 



if the amplitude A in Equ. (13) is chosen according to 



A = ^(^ 



1/4 



(15) 



where is the size of the system in x— direction and 
Rq is defined in appendix A. The gap, for the portion of 
the Hc2-cuvve characterized by n, is a linear combination 
of all (f)n,k belonging to this n. The harmonic oscillator 
cigcnfunctions (13) arc extensively used in the theory of 
the quantum HalP^ effect and many other topics in the 
quantum theory of a charged particle in a magnetic field. 



III. FREE ENERGY EXPANSION NEAR THE 
UPPER CRITICAL FIELD 

We assume that the transition between the supercon- 
ducting and normal-conducting states at the upper crit- 
ical field Hc2 will be of second order for arbitrary tilt- 
angle 9. Then, the order parameter A, or more precisely 
its amplitude e may be used as a small parameter for ex- 
panding the free energy G in the vicinity of Hc2 ■ We keep 
terms up to fourth order in e and all orders in order pa- 
rameter derivatives and determine the energetically most 
favorable order parameter structure near Hc2- Similar 
calculations for the ordinary vortex lattice, correspond- 
ing to the case of large 6 of the present arrangement, 
have been performed by Eilcnbcrger^'^ and by Rammer 
and Pesch-^'*. No special assumptions on the order pa- 
rameter structure, such as the number of zeros per unit 
cell, will be made. We only assume that the order param- 
eter is quasi-periodic on a 2D lattice, with an arbitrary 
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unit cell, characterized by the length of the two basis vec- 
tors and the angle between them. The free energy will 
be minimized with respect to these unit cell parameters. 

A. Order parameter 

Let the unit vector a of our elementary cell be paral- 
lel to the X— axis, a = acx- The angle between a and 
the second unit vector b is denoted by a. To construct 
a quasi-periodic order parameter near Hc2, exactly the 
same method as used by Abrikosov^''', for the case n — Q, 
may be applied. The result is given by the following lin- 
ear combination of a subset of the basis functions (13) 

^n(f) = exp I — t7r-m(m + 1) cosa 1 x 

m= — oo ^ ^ 

exp ^i— mx^ /i„ (y — m6sina) , (16) 

where 

K{z) = I^e-^^'ife„ (^251^) . 
Vn! ^ ' 

This order parameter^'^^'^^ is not invariant under trans- 
lations r ^ r = r + na + mh but acquires phase factors 
for each elementary translation, which are uniquely de- 
fined within a fixed gauge. Surrounding a unit cell in 
anti-clockwisc direction, these phase changes add up to 
a total factor of expi27r, i.e. each unit cell carries a sin- 
gle flux quantum $o • We shall use this assumption of 
a single flux quantum per unit cell, which is written as 
B ^ah'svua = 27r in the present imits, throughout this 
paper. Preliminary calculations^^ show that states with 
two flux quanta per unit cell have higher free energy and 
can be excluded. Also, a preference for multi-quanta vor- 
tices seems unlikely in the present situation, where the 
single flux quantum state is stable at large 0, while the 
total flux decreases to zero as 9 — > 0. 

The order parameter (16) describes a flux line lattice 
where the Cooper pair states belong to arbitrary Landau 
quantum numbers n, depending on the tilt angle G. As 
is well known, the pairing states for the ordinary vortex 
state belong to the lowest Landau level n = 0. The 
present shift to higher Landau levels is, of course, related 
to the large paramagnetic pair-breaking field as will 
be discussed in more detail in section VI. 

The coefficient C„ in Eq. (16) may be expressed by 
the spatial average of the square of the order parameter, 
using the relation 

-FM/2 

(|An|2) = ^ 1, (17) 

^ m=-M/2 

where is the area of the sample. The spatial aver- 
age over the unit cell area = a6sina is defined in 



appendix A. For later use, when performing the limit 
^ in section VI, we assumed in Eq. (17) that the 
area of the superconducting plane is finite and that the 
number of unit cells in one direction is M. At the end of 
the following calculation, C„ will be fixed according to 
the requirement (|A„p) = 1 and an infinitesimal ampli- 
tude e will be attached in front of each power of A„. 

A useful quantity is the square of the order parameter 
modulus, which may be written in the form 

\U\A = Y.^^l\i^'^''^'- (18) 

The Fourier coefficients {il>n)i,j are given by 

= (-l)'^c-™'^ ^°^"e--'-/2L„(xz,,), (19) 

where xij is defined by Eq. (B4). The order parameter 
tpn is proportional to A„ but with an amplitude chosen 
according to (IV'nP) = 1- It is instructive to compare 
Eq. (18) with the local magnetic field reported later in 
subsection III E. 



B. General aspects of the expansion 

A fourth order expansion of G requires first and third 
order contributions in the Greens functions /, /+. We 
use the notation 

f = f^'^+f^'\ r = r^'^+f+^'\ (20) 

where /^^^ and /^^^ are the contributions of order and 
respectively. A consistent treatment of the magnetic 
field terms^*'^^ requires a separation of B and A accord- 
ing to 

B{r) = B + Bi{r), A{r) = I{r} + Ai{r), (21) 

where B is the spatially constant magnetic induction, 

Bi{r) is the r— dependent deviation from B and A(f), 
Ai (f) are the corresponding vector potentials. An evalu- 
ation of the magnetic field terms in G requires the lead- 
ing order in Bi{r), which is e^: i?i « B^'f^ . The spatially 

constant quantity Hc2 — B, where B = \B\,is small of or- 
der e^. The whole expansion in e will be done keeping B 
fixed; at the end of the calculation, the Gibbs free energy 
G will be minimized with respect to the order parameter 
amplitude e and the induction B. The calculation can 
be seen as an extension of Abrikosov's classical work^^ 
to arbitrary temperatures below T^. 

Let us choose the coordinate system in such a way 
that the magnetic field lies in the (y,z)— plane. Then, 
the induction B{r} (and the external field H) may be 
split according to 

B{f) = S|| {f)ey + B^{r)e,, (22) 
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in perpendicular and parallel components B±, The 

corresponding vector potentials are denoted by A±, 
In order to fix the gauge we may employ here essentially 
the same method as used before in numerical calculations 
on the vortex lattice without Zeeman coupling^^'"^". The 
gauge conditions which fix are given by^* 



df 



0, / (frA\ =0, Ai periodic. 



(23) 



The vector potential A describing the average value B of 
the induction is chosen according to 



(24) 



The first term in Eq. (24) can be omitted in the gauge 
invariant derivatives of Eq. (2) since no .2— dependence 
exists in our 2D system. Thus, the orbital pair-breaking 
contribution in the transport equations consists of the 
sum of the second term oc B±_ in Eq. (24) and the 
r— dependent part Ai (only the perpendicular compo- 
nent of A-i is relevant here). The (large) parallel com- 



ponent , on the other hand, enters the spin pair- 
breaking term, which is proportional to B{r) = (i?y (r) + 
B\{r)Yl'^. Eqs. (23,24) fix the gauge, i.e. allow a 
imique determination of A in terms of B. While |Ap 
and B are periodic, i.e. invariant under translations be- 
tween equivalent points in the 2D structure, A and A are 
only quasiperiodic, i.e. they differ by phase factors and 
a change in gauge respectively. The phase factors are 
fixed within a given gauge and may be calculated using 
Eq. (24). 

As a first step in the expansion of G, the Green's fimc- 
tion g is eliminated in favor of /, /+ by means of the 
relation 



5 = 1- 



+ .. 



which is valid for small A. Second, the gradient terms 
in G may be eliminated with the help of the transport 
equations (2). Then, the (dimensionless) Gibbs free en- 
ergy takes the form 



a 



--^ a7+ + av + -(a7(7+F+a*757+)+c.c. 



(25) 



where the bar denotes a Fermi surface average as defined 
in appendix A. 

Inserting the expansions (20), (21) in the free en- 
ergy (25) and collecting terms of the same order in e, 
G takes the form 



^ = (5 + g(2)+gW, 



(26) 



where the terms G, G^^^ , and G^^^ denote the free energy 
contributions of order e*^, and respectively. The 
term G is given by 



2 52 



(27) 



We will first simplify the quantities G^^^ and G^^^ and 
then calculate the minimum of G with respect to the 
amplitude e and the induction B. 



C. Second order contribution 

The second order contribution to the Gibbs free energy 
is given by 

OO 

-i5](A7T(T)+A*7W+c.c.) 



(28) 



To calculate the lowest order Green's function only con- 
tributions of order e°, namely the spatially constant part 

1/2 



B 



( 



B'i + B-i 



of the induction and the lowest or- 



der vector potential A, have to be taken into account in 
Eq. (2). The resulting equation for /^^^ is given by 



[coi + ifiB + kdl^^^jf = A, 



(29) 



where dr = ■§p + iB^yCx. To proceed, we use well- 
known methods"^^ and solve first the eigenvalue problem 



8 



of the operator The solution is given by 

kdi°^fkpi^ = \rh,p(^^ (30) 
with the eigenvalues Ej^ p = ikp and the eigenfunctions 



i=, xy _^ 
zij_L — + ipr 



(31) 



Using the completeness of this continuous set of eigen- 
functions the differential operator on the l.h.s. of Eq. (29) 
may be inverted and /^^^ be represented in the form 



ivi + tfiB + ikp 



(32) 



Representing the denominator in Eq. (32) by means of 
the identity 



from appendix B. The remaining summations and inte- 
grations may be performed analytically^^. Collecting all 
terms one obtains the final result for the second order 
contribution 



=(|A|2) 



Jo 



sinh st 

cos{iiBs)e-^^''/^Ln{B_LS^/2)] 



[1- 



(37) 



While the order parameter expansion Eq. (16), which 
entered the calculation of G^'^\ depends on the lattice 
parameters a, b, a, this dependence is absent in the fi- 
nal result, Eq. (37). The quantity G^^^ characterizing 
the appearance of the superconducting instability, and 
not the detailed structure below it, docs only depend on 
the eigenvalue n. The relation G^^' = agrees with the 
linearized gap equation (12) used to calculate Hc2- 

The technique used here to calculate G^^^ will be gener- 
alized in the next subsection to evaluate the fourth order 
contribution to the free energy. 



dse 



(33) 



as an additional integral, both the p— integration and the 
ri— integration may be performed analytically and the 
solution of Eq. (29) takes the form 



Jo 



f 



exp 



due 



B 



(^-2uykx + u^kxky^ 



A*(r - uk). 



The first order solution for /+ is given by 
/+W(fc,a;„r) = /«*(-fc,<,rl. 

The evaluation of the remaining integrals may be 
greatly simplified by introducing the gap correlation 
function V{fi,f2). In the present gauge it is defined by 



V{n,f2) = A(fi)A*(f2)exp 



B 



r(a;i -X2){yi + y2) 



,(35) 

Of particular importance are the Fourier coefficients 
Vij{r), where r = ri — f2- The precise definition and 
calculation of Vij{r) is reported in appendix (B). 

All terms in Eq. (28) containing first order Green's 
functions may be expressed as integrals over a gap cor- 
relation function. The first of these takes the form 



j-C 

A/+(i) = / 
Jo 



due-'"^'V*{f+uk,r), (36) 



while the corresponding expression for A*/*^^^ may be 
derived from (36) with the help of the symmetry rela- 
tions (7). To proceed, center of mass coordinates are in- 
troduced and a Fourier expansion of y C'*^ (R^ f) with re- 
gard to the variable R is performed, using the result (B3) 



D. Fourth order contribution 

The free energy contribution of order may be split, 
according to 



gw = gW + gW, 



(38) 



in a nonmagnetic part G^'' and a magnetic part G^'' . In 
(34) G^^ the spatially constant induction B and the corre- 
sponding vector potential A{f} are used. The term G^^ 

collects all terms of order e'^ where deviations Bi (r) w 
(or the corresponding vector potential Ai{r)) from the 

average induction B are taken into account; for k — > oo, 
it becomes negigibly small. 



The nonmagnetic part G^^ is given by 



The term G^^ may be calculated using the solutions 
of order e^, already obtained in subsec- 
tion IIIC, 



(39) 



Nd 

E 

1=0 



Ja*/(1)2/-F(1) + A/(l)/+(l)2 + C.C.] ). 

(40) 



The term requires the nonmagnetic parts f^^ , 
of the third order Greens functions f^^\ f~^^^\ 



f+(3) 
JN 



G 



(4) 



1=0 



(3) 
N 



The magnetic part G^^ is given by 



M 

(4) 



Gj^' = GW+G: 



(4) 



(41) 



(42) 
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The term ci^^ is purely magnetic in origin, while the 
term G^'*'' contains the magnetic parts /m'^'' of the 
third order Greens functions, 



(43) 



In a next step, the terms Bi and /'^"'' = /,!} +/m of or- 
der and f"* respectively, must be calculated. The same 
method used in subsection IIIC to calculate /'^', by in- 
verting the differential operator on the l.h.s. of Eq. (29), 
may be used here to obtain /^^^ . Using an operator nota- 
tion for brevity, the sum of f^^ and may be written 



as 



= [u;i + ifiB + kdiP^]-^D, (45) 
D = -Ia/Wz+W-P/W. (46) 



The first and second term in Eq. (46) gives /^'and 
respectively. The term P is of order and is given by 



ikA\ 



•(2) 



(47) 



sj^^-* and A^^-* must be 



The magnetic contributions Bi 
determined by solving Maxwell's equation (5). Expand- 
ing (5) one obtains two decoupled equations 



I, 1 B.„ 



(48) 
(49) 



for the parallel and perpendicular component and 

Bi± of Bi. The quantities 77°, /?, which are both of order 
e^, are given by 



1=0 



i=0 



The Green's fimction g in Eq. (50) may be replaced by 
— /^^'/''''•^V^ under the A;-integral. Thus, 77°, /3 may be 
calculated by using the first order solutions f^^\ f~^^^^ 
as given by Eq. (34). The solution of Eqs. (48, 49) 
is obtained by expanding the unknown variables 

Bi± and the parameters 77°, (3, which are all invariant 
under lattice translations, in Fourier series; the corre- 
sponding Fourier coefficients are denoted by (-Bi||)i.m = 

The explicit solutions will be reported at the end of this 

subsection. 



Given the second order contribution Bi , the first term 
Gc^^ of fl^^ [see Eq. (43)] can be evaluated. To calculate 
the second term G^^^ one needs, in addition, the correc- 
tion term ^1 [see Eqs. (45)-47)]. Writing A 1 = Am+Ai±, 

the Fourier coefficients of Ai± may be expressed^^ 
in terms of the Fourier coefficients of the induction, 



(51) 

(^l±)j,m = -^7— (i?l±);,m {Ql Ql,m,xey) , 



using the gauge conditions defined by Eq. (23). The 
quantities Qi,m,x, Qi,m,y in Eq- (51) are the x and y com- 
ponents of the reciprocal lattice vector Qi^m defined in 
appendix B. 

Each one of the four terms of order e'^ in Eqs. (39, 42) 
may be represented as a multiple integral and Matsubara 
sum over the product of two gap correlation functions. 
What remains to be done is to perform analytically as 
many integrations as possible. The details of the calcula- 
tion will be reported here for the first term gI^' , defined 
by Eq. (40); the evaluation of the other three terms is 
similar. 

Using the first order Green's functions (34) and the 
definition of the gap correlation function (35), the term 

g1^^ takes the form 



J 



fOO pOO fOO = ; 

Gi'^) = -(-V / ds / dsi / ds2e-'^^^'+''+'^'>V{f- sik,f+S2k) V{f- sk,r) + V{f,f+sk) 

° Jo Jo L 

Expanding F in a Fourier series, the spatial average in Eq. (52) may be performed and g'^^ takes the form 
Ga^=-oY.Y E / / dsi / ds2 e-'(''+«^+«^)l^,„ (-(si + s^) Vi*^ (sk) cos (Qi,^-k 



+ c.c. 



)■ 

(52) 



cos ( cj/,m— 'iT-^fc ) cos {iJiB (s -I- Si -I- 82)) + sin ( Qi 



^-fc ) sin (/iS (s -I- si -h S2)) 



where Vi^m is given by Eq. (B3) and the symmetry relations (7) have been used to rearrange the integrand. We 
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introduce center of mass coordinates ts = si + S2, tji = si — S2 in the si,S2-plane. Replacing si,S2 by the new 
variables, the integration over tji may be performed and the double integral over s and ts becomes a product of two 
independent, one-dimensional integrals. Performing this step, G^^ takes the form 

Nr 



^i'^ = - * E ^ / E / da e-'^ cos (-Qi,mk) ■ (sk) . 

7^0^'' Jo i^^Qi,mkJo \2 J \ J ^^^^ 

dts e~"'*® sin ( -^Qi,mk \ ■ Vi^m (^-tskj [cos (fiBs) cos {nBts) - sin (/zSs) sin [fiBts)] ■ 



f 

Jo 



An attempt to simplify Eq. (53) further, by performing 
one of the remaining integrations analytically, was not 
successful. At this point it seems already feasible to cal- 
culate the remaining integrals over s, (p and the sums over 
Matsubara and Fourier indices numerically. However, we 
prefer to proceed and calculate the remaining integrals 
by means of an asymptotic approximation. 

Let us consider for definiteness the integral over 
s in Eq. (53). The integrand has its maximum at 
s = 0. We analyze the behavior of the various fac- 
tors in the integrand as a function of s, and ne- 
glect the s— dependence of the slowest varying factors. 
The characteristic lengths in s— space of the factors 

exp(— w;s), cos^sQi fnkj: Vijn(slt) and cos (/i_B,s) are 
given by n = [{21 + l)t]-\ = iB^\f_i^rn\)~\ Js = 
(n_Bj_)~^^^ and T4 = {^B)~^, where fim = la + mb. We 
consider a range of inductions B ^ Bp, where the Pauli 
critical field Bp is (in the present system of units) given 
by ^Bp = 0.4. As a consequence T4 > 2. Choosing a 
typical number /i = 0.1 for the dimensionless magnetic 
moment, our induction varies in the range B ^ 4. The 
characteristic lengths T2 and T3 both depend on the Lan- 
dau quantum number n: recall that B± depends on n as 
shown in Figure (1). Let us consider first the case n — 1. 
Then, B± = BsinOi with sinGi Ri fx/n according to 
Eq. (1) and the definition of /i in appendix A. As a con- 
sequence, T3 is of the same magnitude as T4 for n = 1. 
The magnitude of of T2 = 7/|n^m| varies strongly depend- 
ing on the Fourier indices l,m. For not too large Fourier 
indices and nearly all uii, t\ will be the smallest of the 
four characteristic length. This is, however, only true for 
not too low temperatures t. For large Fourier indices, 
which should be taken into account in the present situa- 
tion, the behavior of the integrand will be dominated by 

the term cos (^sQi^mkj because its characteristic length 

T2 becomes small for large l,m. Thus, the latter term as 
well as the Matsubara term exp(— wjs) has to be kept, 
while the terms Vi^rn{sk) and cos (/iSs) show the slow- 
est variation in s and may be replaced by their values 



at s = 0. This conclusion remains true for arbitrary n. 

This may be seen by using the relation nB± = /3 which 
will be derived in section VI. 

Using this asymptotic approximation both the integral 
over s and the Fermi surface average may be performed 
analytically and one arrives at the result 

1=0 l,m 

24 + j\Q,.n,\- 

The second nonmagnetic term G^^ [see Eq. (41)], which 
is evaluated with the same method, is given by g'^^ = 

In order to calculate the fourth order terms of mag- 
netic origin, Gc^' and G^f^ , the Fourier coefficients of the 
quantities 77", f3 [see Eq. (50)] have to be evaluated first. 
This may be done using a method similar to the one out- 
lined above for G'^K A noticeable difference is, that the 
s— dependence of the slowest varying factors (the ones 
with characteristic lengths T3 and T4) cannot be com- 
pletely neglected in the course of the asymptotic approx- 
imation, but must be taken into account to linear order in 
s. In a second step. Maxwell's equation has to be solved 
to obtain the magnetic field correction Bi. Given the 
latter, the free energy Gc^^ may be calculated. The term 
G^'*'* contains an additional s— integral which may be per- 
formed by means of an asymptotic approximation of the 
above type. The relation G^^^ = —2G'^^ was again found 
to be true, in analogy to the nonmagnetic case (the same 
relation has been found in microscopic calculations^^ of 
the ordinary vortex lattice near Hc2)- 

Collecting all fourth order terms and attaching the fac- 
tor one obtains the final result 



G(4) 




(1) 



l.m J 
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where the prime at the second summation sign indicates 
that the term ? = 0, m = is to be excluded from the 
sum. The functions /i , gi depend expUcitly on the Lan- 
dau quantum number n and are given by 



fi{x) = e ^Ln{x) 



e 2 



l-Lnix) + {1 - Sn,o) Ll_,{x) 



(55) 
(56) 



where L\ is a Laguerre polynomiaP^. The Matsubara 
sums are given by 



Nd 



-AQi 



1=0 ojf (wf + \\Qi,m\^) 



3/2 • 



o(2) _ V- 



1/2- 



(57) 

(58) 



The square of the reciprocal lattice vector is conveniently 
written in the form |(5i,mP = 2i?_LXz_m/2 where xi^m is 
defined by Eq. (B4). Introducing a magnetic length L 
defined by 



B± = 



27r 



2 

Z2' 



(59) 



ab sin a 

these parameters which depend on a, b, a, are given by 



Sin a \a 



L 



9 „ , cos a 

-) TO^-27r/m . (60) 

\L/ sma 



E. Local induction 

The components Bi± and of the spatially varying 
magnetic field Bi are given by Fourier series of the form 



SiA(r)=X^(BiA)z,me*'5'-^ 



(61) 



l,m 



where A =_L, ||. The Fourier coefficients are given by 



_ t(|A„P) , 

{Bl\\}l,m - - ^2 72^\\t^ ■ 

r\t jji 

(-1) e - fl(Xl,m)bi^^, 
(i<l±)i,m - - -2-(-l) e 

rv fJi 



(62) 



{B^t^'fi{xi,m)Slll - 9i{xi,m)SlX) . (63) 

The parallel component (62) is proportional to /it^ and 
is entirely due to the spin pair-breaking effect. The per- 
pendicular component (63) is the sum of a /i^-dependent 
term and a second term not (explicitly) dependent on /x. 
The terms dependent on /i^ have the same form for both 
components (recall that the direction of B in spin space 



is arbitrary) and arc proportional to the relevant compo- 
nent of the macroscopic induction. The second term in 
Eq. (63), which is of opposite sign, may only forn = be 
considered as a consequence of orbital pair-breaking; for 
n > this second term depends also (since a positive n 
is necessarily due to a finite /i) on the spin pair-breaking 
effect. The GL limit of the local induction is discussed 
in appendix C. 

The validity of the asymptotic approximation used in 
the derivation of Eqs. (54,61) is not restricted to low n, 
but sufficiently high temperatures, say t > 0.1, should be 
used. Clearly, if different states with very small free en- 
ergy differences are found, no conclusion as to the relative 
stability of these states can be drawn. 



F. Extremal conditions 

In thermodynamic equilibrium, the values of e, B±, By 
and the lattice parameters a, 6, a have to be chosen in 
such a way that the free energy becomes minimal. To find 
the equilibrium values of e, B±, B|| the extremal condi- 
tions 



dG 



= 0, 



dG 

mi 



= 0, 



dG 
dBl 



= 0, 



(64) 



have to be solved near Hc2 ■ The question for the optimal 
a, b, a will be addressed in the next section. 

Inserting the superconducting solution for e in the free 
energy yields 



G = G 



1 (G(^)) 
4 G(4) 



(65) 



where the coefficients G G"^, are defined by G = 
G + e^G^^) 4- £4^(4). Eq. (65) shows, that the stable 
lattice structure (see section IV) is determined by the 
requirement of minimal G4 

To find the two-component macroscopic magnetiza- 
tion relation between induction B_\_, B|| and external field 
i?_L,i?l|, the above extremal conditions must be solved 
for B±, B\\. This cannot be done for arbitrary fields 
but requires an appropriate expansion of the coefficients 
for small B±_ — Bc2,±, -B|| — Bc2,\\- A lengthy but 
straightforward calculation, generalizing Abrikosovs clas- 
sical work^^ to the present situation, leads to the result 



B± = aj_±Hj_ + Q;_L||-ff|| + P± 
=a±l|i?± + aim 4-/3,1 . 



(66) 
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The coefficients in this Unear relation are given by 

aj__i_ = 2k^ [2 (k^ - n'^) - A\\\ /detM 

o^^\ = 2k'^A\\^/detM 

ami = [2 _ _ /detM 

/3± = -2 (k' - m') (^±Se2,± + ^||±S,2,||) /detM 

/3|l = -2 (k^ _ (Aj|Be2,ii + A\\^Bc2.,i) /detM, 

where: detM = 2 (k^ - y^^) [2 (k^ _ fj^) _ _ . 

The parameters ^y,... may be calculated for a given 
lattice structure with the help of the relations 



An = 



1 



2G(4) ^ dB\\ 



1 ( d&^ 



1 9G(2) ^^(2) 
2(5(4) dBu dB_L 



(67) 

(68) 
(69) 



where the derivatives of G^^) have to be evaluated at 
B = B(.2 and the relation = ^yj^ may be shown 

to be true. 

Eq. (66) constitutes the macroscopic relation between 
induction and external field for a 2D superconductor in a 
tilted magnetic field. It is, of course, strongly anisotropic 
and shows a coupling between the parallel and perpen- 
dicular field components. For = 0, and t ^ 1, 
Eq. (66) should reduce to Abrikosov's GL solution^^'^^ 
for the magnetization of a triangular vortex lattice. This 
is indeed the case as shown in Appendix C. 

For = 0, /U 0, Eq. (66) describes the ordi- 
nary vortex lattice (near Hc2) for arbitrary tempera- 
tures. A numerical comparison with corresponding re- 
sults by Eilenberger^^ and Rammer and Pcjsch^^ has not 
been undertaken because a difi'erent (spherical) Fermi 
surface has been used in these works. However, the limit 
= 0, /X of the present theory will be checked in 
appendix D by calculating the critical value of k separat- 
ing type II from type I superconductivity. 



IV. RESULTS FOR FINITE PERPENDICULAR 
FIELD 



In this section we determine the stable order param- 
eter structures for the paramagnetic vortex states with 
1 < n < 4 in the vicinity of Hc2- The numerical pro- 
cedure to find the stable states is essentially the same 
as in KRS®. First, the upper critical field Bc2 and the 
corresponding quantum number n have to be found for 
given temperature t and tilt angle O by solving the lin- 
earized gap equation (37). In a second step, the stable 
lattice structure, which minimizes the fourth order term 
^(4) ^ G(4)/g4 jgg^ (54]^ Yias to be determined. Be- 
cause of the fiux quantization condition the minimum 



with respect to only two parameters, which may be cho- 
sen as a/L and a, must be found. In contrast to the 
ordinary vortex lattice, where it is usually sufficient to 
calculate only a few lattices of high symmetry (triangu- 
lar, quadratic) to find the stable state, the present situ- 
ation is characterized by a large number of local minima 
of Eq. (54), corresponding to a large number of possible 
lattices of rather irregular shape. Therefore, a graphi- 
cal method was used to determine the stable state; the 
free energy surface {a/L, a) was plotted for the whole 
{a/L,a)-pla.ne and the global minimum was determined 
by inspection. Basically, two material parameters, fi 
and K, and two externally controlled parameters, t and 
G, enter the theory. Numerical calculations have been 
performed for a single value of jj. = 0.1, two different 
reduced temperatures 0.2 and 0.5, four different values 
0.1, 1.0, 10, 100 of Eilenbergers parameter K, and several 
values of 6 corresponding to different Landau quantum 
numbers n. Some of the resulting order parameter and 
magnetic field structures in the range n < 4 will be re- 
ported here. These low-n pairing states are, of course, 
the most important ones from an experimental point of 
view. 

For comparison we consider first, in appendix D, the 
ordinary vortex lattice state with n = 0. This illustrates 
the method and may also be used to check the acciiracy 
of our asymptotic approximation. The equilibrium state 
for \ow-K type II superconductors is calculated and good 
agreement with previous theories is found for not too low 
temperatures. 

Considering now pairing states with n > 0, the num- 
ber of order parameter zeros per unit cell increases 
clearly with increasing n. One finds* two types of min- 
ima of &^\a/L,a), isolated minima and line-like min- 
ima. The first type corresponds to "ordinary" 2D lat- 
tices, the second type, characterized in a contour plot 
[see Fig.l of KRS*] by a line of constant a/L with 
&'^\a/L,a) nearly independent of a, corresponds to 
quasi-one-dimensional, or "FFLO-like" lattices (rows of 
vortices and one-dimensional FFLO-like minima alter- 
nating). A convenient way to identify the type of 
minimum and find its position on the a/L— axis, is to 
plot the projection of the G^^-* (a/L, a)— surface on the 
(G^^-*, a/L)— plane. An example for this perspective, 
where a— independent parts of the free energy surface 
show up as lines, is given in Fig 2 for n = 7. The 
a— coordinate of a 2D minimum cannot be read off from 
such a plot and requires a second projection on the 
(a/L, a)-plane (such as Fig 11 or Fig.l of KRS*). The 
free energy maps for other n > states are in principle 
similar to Fig 2 but the different local minima show more 
pronounced differences for smaller n . 

Let us start with the paramagnetic vortex state with 
n=l and consider first the limit of large n. As reported in 
KRS*, a quasi-one-dimensional state is found to be sta- 
ble in this case. Fig 3 shows the spatial variation of the 
modulus of the order parameter. One sees rows of vor- 
tices separated by a single, FFLO-like line of vanishing 
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I ' 1 ' 1 ' 1 ' 1 < 1 — 

0.40 0.80 1.20 1.60 2.00 2.40 

a/L 



FIG. 2; (Reduced image quality due to arXiv restrictions) 
Projection of the free energy G^*' on the G(*\ a/L-plane. 
Using this perspective the a— independent parts of the free en- 
ergy surface are displayed as lines. In the considered range of 
a/L one finds six local minima, corresponding to two FFLO- 
like and four two-dimensional lattices. The global minimum 
is at a/L ~ 1.1 and corresponds to a two-dimensional lat- 
tice. Parameters chosen in this plot axe n = 7, t = 0.5, R, = 
10, fi = 0.1, e = 0.055. 



order parameter. The unit cell of the structure shown in 
Fig 3 is given by a/L = 1.0875, a = 33deg. A shift of 
the vortex rows relative to each other leads to a lattice 
with the same a/L and a different a, which has nearly 
the same free energy (which is reasonable, since the inter- 
action between vortices from different rows is weak as a 
consequence of the intervening FFLO domain wall). The 
vortices are of the "ordinary" type, i.e. the phase of the 
order parameter changes by +2^ when surrounding the 
center. 

It is of interest, to calculate the magnetic field belong- 
ing to this order parameter structure. We plot the par- 
allel and perpendicular components i3i||(r) and Bi±{r) 
of the spatial varying part Bi{r) of the magnetic field 
as given by Eqs. (61)-(63), omitting a common factor 
i(|A„|2)/(K2 - /i2). The field Bi||(f), which is entirely 
due to the spin pair-breaking mechanism, is shown in 
Fig 4. Due to its paramagnetic nature, the field Bi\\{r) 
is expelled from regions of small This behavior is 

exactly opposite to the usual orbital response, which im- 
plies an enhancement of the induction in regions of small 
I'f/'Kr). As a consequence, the spatial variation of is 
very similar to that of ji/'P, shown in Fig 3. 

The perpendicular field Bi±(f), shown in Fig 5, con- 
sists of a spin term proportional to /i^, and a second term 
which depends [see Eq. (63)] not explicitly on fi. The 




FIG. 3; (Reduced image quality due to arXiv restrictions) 
Square of modulus of order parameter as a function of 

x/a, y/a in the range < x/a < 2, < y/a < 3.2. This is the 
stable structure (unit cell parameters a/b = 0.205, a = 33 ° ) 
for t = 0.2, k = 100, e ^1.2° (n = 1, Bc2 = 4.141) 




FIG. 4: (Reduced image quality due to arXiv restrictions) 
Parallel component Biy as a function of x/a, y/a in the range 
< x/a < 2, < y/a < 3.2. This plot has been produced 
using the same input parameters as in Fig 3. 



14 



term proportional to fM^ is negligibly small and the total 
field is essentially given by the second term. Near the 
vortices the field Bi±{r) behaves in the familiar, orbital 
way, i.e. it is largest at the points of vanishing ijj and de- 
creases with increasing distance from the vortex centers. 
However, at the FFLO-like lines of vanishing order pa- 
rameter, where no topological singularity occurs, Bi± has 
a minimum, i.e. shows paramagnetic behavior. Thus, the 
magnetic response of a n = 1 superconductor may either 
lead to a local suppression or to an enhancement of the 
magnetic field in regions of small order parameter. This 
is in contrast to the purely orbital response of a n = 
superconductor, where the magnetic field is always en- 
hanced. This unconventional behavior is formally due to 
the second term in gi [see Eq. (55)]. The field is 
much smaller than Bi± and the total field Bi for n — 1 
is consequently dominated by the perpendicular compo- 
nent Bi±, which is a consequence of the combined action 
of both pair-breaking mechanisms. 

The quasi-one-dimensional order parameter structure 
shown in Fig 3 seems to be representative for the pairing 
state with n — 1; no other stable state has been found for 
K = 10 and t — 0.5. At k = 1, 0.1 the free energy surface 
has no minimum at all, which means that a transition 
to type I superconductivity occurs at some value of k 
between 1 and 10. 

Extrapolating the n — 1 result to higher n, one would 
expect the following structure for the pairing state with 
Landau quantum number n: rows of vortices separated 
by n lines of vanishing order parameter . Such a struc- 
ture would approach the (line-like) FFLO state in the 
limit n =4> oo. However, this simple picture is not real- 
ized, at least in the important range of low n. It holds 
generally for odd n, but for even n two-dimensional struc- 
tures are preferred. In the latter case, one has n + 1 
isolated order parameter zeros per unit cell, with associ- 
ated phase changes of a multiple of 27r. Such a situation 
leads necessarily to the presence of one or more antivor- 
tices - vortices with a topological phase change of —2tt 
around the center - for states with even n, since the to- 
tal phase change around the unit cell must remain +2tt. 
Recently, various proposals to create stable antivortices 
have been published; see e.g. Moshkalkov et al'^^. In 
the present context, it is clearly the strong paramagnetic 
pair-breaking, which is responsible for the stability of the 
antivortices. Among the (even-n) antivortex states, the 
one with n = 2 is most easily accessible from an experi- 
mental point of view and very stable under variations of 
t and K. Its properties will be discussed in detail in a 
separate publication'^^; a preliminary account has been 
published already'^^. 

For n = 3 free energy minima for spatially varying 
states exist in the whole considered range 0.1 < k < 
100 of the GL parameter. Thus, increased spin pair- 
breaking stabilizes inhomogeneous equilibrium structures 
and shifts the phase boundary between type II and type 
I superconductivity to lower values of k. In the high-K 
region (for k > 10) the stable state of a n = 3 supercon- 




FIG. 5: (Reduced image quality due to arXiv restrictions) 
Perpendicular component Bi± as a function of x/a, y/a in 
the range < x/a < 2, < y/a < 3.2. The same input 
parameters as in Fig 3 have been used. 



ductor is of the quasi-one-dimensional type (at lower n a 
2D state of nearly the same free energy has been found, 
which will not be discussed here). The fields -Bi|| 
look similar to the n — 1 case (see Figs. 3 and 4) except 
that the vortex rows are now separated by three FFLO- 
like lines of vanishing order parameter. The vortices in 
neighboring rows are already completely decoupled for 
n = 3; a translation of neighboring rows relative to each 
other changes the angle a between the unit cell basis 
vectors but does not lead to any change (within 8 dig- 
its) of the free energy. The perpendicular induction Bi± 
is again dominated by the second term in Eq. (63) and 
looks similar to the n — 1 case (see Fig. 5); in contrast 
to the spin part this field does not reflect the detailed 
order parameter structure but has only a single broad 
minimum at the position of the three FFLO lines. 

For a n — A superconductor a,t t — 0.5 the equilibrium 
state is of the quasi-one-dimensional type for k < 1, and 
of the 2D type for k, > 10. Fig. 6 shows the 2D order 
parameter structure for a superconductor with k = 10. 
There are 5 zeros of ^p per unit cell, one of them of 
an elongated shape. The nature of these topologically 
singular points may be clarified by plotting either the 
phase'^^ or the local magnetic field. The parallel compo- 
nent -Bi||(r) of the field i?i(r) is again (compare Figs. 3 
and 4) similar in shape to the order parameter \tp\'^ and 
need not be displayed here. The perpendicular field Bi± 
is shown in Fig. 7. Three of the 5 order parameter zeros 
displayed in Fig 6 belong to "ordinary" vortices, with lo- 



15 




FIG. 6; (Reduced image quality due to arXiv restrictions) 
Square of modulus of order parameter \tp4\^ as a function 
of x/a, y/a in the range < x/a < 1.4, —0.2 < y/a < 
0.62. This is the stable structure (unit cell parameters a/b — 
0.6735, a = 70.125° ) for t = 0.5, k = W, 9 = 0.1 ° (n = 
4, Bc2 = 3.486) 




FIG. 7: (Reduced image quality due to arXiv restrictions) 
Perpendicular component Bi± as a function of x/a, y/a in 
the range < x/a < 1.4, -0.2 < y/a < 0.62. The same 
input parameters as in Fig 6 have been used. 



cal field enhancement and diamagnetic screening current 
(two of the three maxima of Bij_ are pronounced, while 
the third, the one corresponding to the elongated zero of 
■0, is rather flat). The remaining 2 order parameter zeros 
belong to antivorticcs with opposite sign of the "screen- 
ing currents" (which are now paramagnetic in nature) 
and with minima of Bi± at the points of vanishing tp. 

Results for n > 4 will not be reported here. Many 
interesting and complex structures may be produced for 
larger n. However, the number of different states with 
similar free energies increases with increasing n. As a 
consequence, the approximate nature of our analytical 
calculation does not allow an identification of the stable 
state for large n . At the same time, an experimental ver- 
ification of these large-n states seems difficult since a very 
precise definition of the tilt angle would be required. 



V. STRUCTURE OF THE FFLO STATE 

The stable state in the purely paramagnetic limit n — > 
oo has been determined first by Larkin and Ovchinnikov^ 
at T = for a spherical Fermi surface. They predicted 
a one-dimensional periodic order parameter structure of 
the form A(r) w cos{qr), which will be referred to as 
LO state. Later, various analytical investigations of the 
stable states in the vicinity of Ttri and near T — 
have been performed'^^''^^ ; many other references may be 
found in a recent review article'^®. A careful search for 
the state of lowest free energy, comparing several pos- 
sible lattices in the whole temperature range, has been 
reported by Shimahara^^. He found that below t = 0.24 
various 2D periodic states have lower free energy than 
the one-dimensional cos{qr} state. Shimahara uses the 
same cylindrical Fermi surface as we do and his results 
do therefore apply to the present problem. Nevertheless, 
we reconsider in this section the problem of the determi- 
nation of the FFLO structure, in order to have a com- 
plete description of all states in a tilted field in a single 
theoretical (quasiclassical) framework. 

The results derived in section 111 cannot be used to 
perform the limit n oo and determine the stable state 
in the purely paramagnetic limit. However, the general 
formalism may be applied in a straightforward way to the 
simpler case of vanishing vector potential. In order to 
be able to compare with previously published results we 
neglect in this section the possibility of spatial variations 
of B and restrict ourselves to the high-K limit. 

The space of basis functions, which has to be used to 
expand all variables near HppLo{T), is now given by the 
infinite set exp(«gr^ with a fixed value of Usually, one 
assumes that the order parameter A fulfills some further 
symmetry (or simplicity) requirements, which then leads 
to a strong decrease of the number of unknown coeffi- 
cients. Following this convention, we restrict ourselves to 
two- and one-dimensional periodic structures. For > 0, 
the order parameter is not periodic but changes its phase 
by certain factors under translations between equivalent 
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lattice points. These phase factors are proportional to 
the perpendicular induction [cf. Eq. (35)] and vanish for 
0^0. Thus, the assumption of a periodic order param- 
eter for = is reasonable (though not stringent). It 
implies, that all allowed wave vectors in the expansion of 
A must be vectors of a reciprocal lattice. 

A further slight simplification stems from the behavior 
of the quasiclassical equations under the transformation 
r => — r, => —k, which implies that the order parameter 
must be either even or odd under a space inversion f => 
—f. Thus, the order parameter may be written as an 
infinite sum 



(70) 



The second order term is given by 

I 



(74) 



with the i— independent coefficient A defined by 
1-e- 



A = 2{\nt + t j 



As- 



sinh st 



■ [l — cos{iJ,Bs)Jo{sq) 



The condition A = determines the upper critical field; 
it may also be derived from Eq. (37), performing the limit 
n — > oo. 

The fourth order term is given by 



with coefficients defined by 



lA 



(71) 



Here, a shorthand notation m is used for the two inte- 
gers characterizing a 2D reciprocal lattice vector Qm [cf. 
the Fourier expansion at the beginning of appendix B]. 
The vectors actually entering the expansion are distin- 
guished by an index i, their total number is /, and the 
two integers characterizing Q„. are denoted by rij. The 
complex numbers Cj are the expansion coefficients; one 
may set ci = 1 since only the relative weight is impor- 
tant. It turns out, that the two solutions distinguished in 
Eq. (71) by a sign are essentially equivalent, and only one 
of them, say the even one, need be considered. Thus, the 
order parameter becomes a linear combination of cosine 
functions. 

All reciprocal lattice vectors used in Eq. (70) must be 
of the same length. Denoting this length by q{T), the 
condition \Qm\ = <l{T) takes the form 



1 



sin Ur, 



1 



Ij cos Up _ 



sin^ Up db sin^ Up 



(72) 



where the two integers I, j have been used here to rep- 
resent the double index m. The dimensionless quantities 
a, b are defined by a = g(T)ap/27r, b = q{T)bp/2Tr, where 
Up, bp, ttp denote the lattice parameters in the paramag- 
netic limit. If / reciprocal lattice vectors exist, the lattice 
parameters Op, bp, ap, fulfill I relations like Eq. (72) with 
I pairs of integers li,ji,... 

Using Eq. (70) the free energy expansion near 
Hfflo{T), including terms of fourth order in the small 
amplitude |A|, may be performed by means of methods 
similar to section III. The result for the purely param- 
agnetic free energy Gp takes the form 



Gp = Gp + G(2) + G^p 



(4) 



(73) 



where G^ 



(2) ^ 

^^H^ , and Gp and G'^^ are contributions 



(4) 



of order |Ap and |A|'' respectively. 



= |Ar[5^ \ci\''A + Y, \ci?\ck?B,^k+ 

I 



(75) 



G^^ depends on the lattice structure via the coefficients 
Ai, Bi^k, and Ci^k, which are defined by 

~ / ^ I T, ( (^) "I" ni (^) ) 



B 



i,k 



Xk) + P, 



ni,—nk,-nk 



Nd „2tt 



(i)) 



where 



Pni,ni,n3 (^) 



U)l +1 



±HB + Qnk 



In contrast to Eq. (54) no approximations have been used 
in deriving Eq. (75). 

Using Eq. (72) all possible 2D lattices and wave vec- 
tors may be calculated numerically. The stable lattice 
at Hfflo{T) is then determined from the condition of 

lowest g'^p^, taking also the LO state into considera- 
tion. It turns out, that it is energetically favorable at 
Hfflo{T) if all eigenfunctions in the order parameter 
expansion (70) have equal weight, i.e. c, = 1 for all i. 

The result of the numerical search for the lowest free 
energy of periodic structures, characterized by maximal 
three pairs of reciprocal wave vectors, is displayed in 
Fig 8. The highest curve at a given temperature cor- 
responds to the stable lattice. For 0.22 < t < 0.56 the 
one-dimensional LO state is realized. For t < 0.22 2D 
periodic structures appear, namely the square state for 
0.05 < t < 0.22, and the hexagonal state for t < 0.05 (we 
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use here the notation of Shimahara'^^ for the 2D states). 
Besides the fact that the triangular state"^^ is absent, 
because it is neither even nor odd, the present results 
agree quantitatively with those of Shimahara^^, obtained 
within a different, but equivalent, formalism. Thus, more 
complicated 2D periodic structures than those found al- 
ready in Ref^^ do not exist in the considered range of 
temperatures; the assumption of equal weight for differ- 
ent wave vectors [cj = 1 for all i in Eq. (71)] has also 
been confirmed for these states. 
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FIG. 8: (Rxxhiccd image quality due to arXiv restrictions) 
The fourth order term Gp*' (minimized with respect to |A|) 
divided by -A"^ [see Eq. (74)] at Hfflo{T) for three differ- 
ent periodic structures as a function of reduced temperature 
t = T/Tc- The part of the hexagonal curve which is lower 
than the LO state is not visible, since the coefficients a are 
determined automatically to yield the highest possible solu- 
tion for -Gi^^/A''. 



The temperature region below t = 0.01 has been inves- 
tigated recently by Mora and Combescot^^. They found 
a series of states characterized by an even (total) num- 
ber 2N = 8, 10, . . . of different wave vectors, all entering 
the order parameter expansion with equal weight, and 
with N increasing with decreasing temperature. Merg- 
ing these results with the present ones, one obtains a very 
simple description of all of the FFLO states at the phase 
boundary, namely an inhnitc number of states, each one 
being a linear combination of iV = 1,2,... cosine func- 
tions of equal weight and with A'' different, but equally 
spaced, wave vectors. 

Of course, it is also of interest to investigate the possi- 
ble equilibrium structures in the region helow the critical 
field. As a first step in this direction, preliminary cal- 
culations at O-'dbMpFLO and Q-'dQEpFLO have been per- 
formed, using the fourth order expansion (73), which is 
not valid near first order transition lines. The result is 
surprising and shows a revival of the LO state in the low 



temperature region, 
0.60 




FIG. 9: (Reduced image quality due to arXiv restrictions) 
The term —G^p^/A^ for the LO state and the square state at 
Q.^QHfflo{T), as a function of reduced temperature T/Tc. 
The hexagonal curve is lower than the LO state and is not 
displayed in this figure. 



In Fig. 9 the terms — Gp'^-'/A^, for the three states dis- 
played in Fig 8, are plotted as a function of temperature 
helow the transition line, at 0.9Hfflo- The hexagonal 
state (not visible) does not exist any more. The usual 
square state (characterized by = 1) is only stable in 
a very small temperature interval 0.061 < t < 0.075. 
The LO state is now stable in a much larger interval 
0.075 < t < 0.56, as compared to Fig 8. It is also stable 
in a small temperature region below t = 0.061. But at 
t w 0.017 the factor -G^^ /A"^ for the LO state has a sin- 
gularity and jumps from -|-oo to —oo. This implies that 
the fourth order term (for the LO state) changes sign 
and that a first order transition occurs somewhere in the 
vicinity of this singularity; higher order terms in the free 
energy would be required for a quantitative treatment. 
Between this singularity at f ~ 0.017 and the lowest con- 
sidered temperature t = 0.01 the stable state is again 
characterized by a square unit cell. However, the order 
parameter in this temperature range, 0.01 < t < 0.017, 
is given by a linear combination of plane wave states [see 
Eqs. (70), (71)] with a real coefficient ci = 1 and an imag- 
inary coefficient C2 = ^. The usual order parameter struc- 
tiirc for the square lattice, which is characterized by two 
real weight factors of equal magnitude (ci = C2 = 1), is 
not equivalent to this case and has higher free energy. 

The results below Hfflo{T) indicate, that the 2D 
states are only stable in a tiny interval near the phase 
boundary, and that the one-dimensional LO state reap- 
pears inside the superconducting state. The square state 
- the one with the smallest N [N = 2) - has the largest 
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stability region, as one would also expect from the free 
energy balance shown in Fig 8. We shall come back to the 
question of the stability of the 2D states in section VI, 
considering it from a different point of view. The struc- 
ture found below the singular point of the LO state (see 
Fig. 9) raises the question, if still other order parameter 
structures, different from those found at the transition 
line, will appear near t = deep in the superconducting 
state. The present fourth order expansion is not really 
appropriate to answer this question. 



VI. TRANSITION TO THE PURELY 
PARAMAGNETIC REGIME 

The limit n ^ oo of the series of paramagnetic vortex 
states, discussed in section IV, is now well known; for 
0.22 < t < 0.56 the one-dimensional LO state is real- 
ized, while 2D states of square or hexagonal type, pre- 
dicted by Shimahara^^, appear at lower t. The region of 
still smaller t, below t = 0.01, which has been studied 
by Mora and Combescot''^, will not be considered here. 
The way, this limit is approached, is, however, unknown. 
Thus, we address ourselves in this section to the the ques- 
tion of how the one- or two-dimensional unit cell of the 
FFLO state develops from the unit cell of the paramag- 
netic vortex states if the Landau level index n tends to 
infinity. 

This limiting process is very interesting, because a 
vast number of different states with different symmetry 
is passed through in a small interval of tilt angles 6. The 
unit cell of the finite-n states is subject to the condition 
that it carries exactly a single quantum of flux of the 
perpendicular field B±. Since B± — » as n ^ oo, at 
least one of the unit cell vectors must approach infinite 
length - i.e. the dimension of the macroscopic sample - in 
this limit. Thus, the n oa limiting process describes 
a transition from a microscopic (or mesoscopic) length 
scale to a macroscopic length scale. 

The transition to the FFLO state has previously been 
investigated by Shimahara and Rainer^ in the linear 
regime. They found the important relation 



q = lim ^/4eB±n/hc, 



(76) 



where q is the absolute value of the FFLO wave vector 
(here we changed to ordinary units). Eq. (76) has been 
derived by identifying the asymptotic form of the Hermite 
polynomials^^ with the form of the LO order parameter. 
It implies that a relation 



n 



a-- 



hcq^ 



(77) 



holds at large n. The validity of Eq. (76) may also be 
checked numerically by comparing the numbers /? and 
q, which are both obtained from the upper critical field 
equation. 



Relation (77) may be derived from basic physical prop- 
erties of the present system. The energy spectrum for 
planar Cooper pairs in a perpendicular magnetic field 
B±_ is the same as for electrons and is given by 

En = nu{n +1), w = — . (78) 
2 mc 

Considering now the energy spectrum of Cooper pairs 
for B± = 0, one has to distinguish two cases. First, in 

the common situation without a large spin-pair-broaking 
field, all Cooper pairs occupy the lowest possible energy 
E = 0, which is the kinetic energy /Am taken at the 
Cooper pair momentum p = Q. Second, if a large spin- 
pair-breaking field parallel to the conducting plane exists, 
the energy value to be occupied by the Cooper pairs, 
shifts to a finite value p^/4m, since the Cooper pairs 
acquire a finite momentum p due to the Fermi level shift 
discussed in section I. Thus, in the latter case, which is 
of interest here, the Landau levels (78) must obey the 
condition 



p 



„ Tie ^ , 1, 
En = —B^{n+-) — 

mc 2 4m 



±r (79) 



for B±_ -^Q. lip is replaced by the wave number q = p/U, 
Eq. (77) becomes equivalent to Eq. (79). The limiting 
behavior expressed by Eq. (76) or Eq. (77) is therefore 
a direct consequence of Landau's result for the energy 
eigenvalues of a charged particle in a magnetic field. 

Combining Eq. (77) with analytical results at T = 0, 
the limiting behavior of the unit cell as n ^ oo may be 
understood. Expressing the FFLO wave number in terms 
of the BCS coherence length by means of the relation 
q = (2/7r)^(^^, and using the flux quantization condition 
in the form 



(80) 



[with $0 = hc/2e and B±^ defined by Eq. (77)] the area 
F^") of the unit cell for pairing in Landau level n is ap- 
proximately given by 



(81) 



Thus, the unit cell area diverges with the flrst power of 
n. The behavior of the magnetic length L , which is 
defined by the relation B± = {^o/'jt)L~'^ , is given by 

Eq. (81) is not sufl&cient to determine the shape of the 
unit cell in the limit of large n. However, a simple possi- 
bility to produce a one- dimensional periodic LO structure 
for n — > 00 is a divergence of one of the unit cell lengths, 
say b, of the form h ^ while the second length a re- 
mains constant, i.e. a « n°. The numerical results for 
the states referred to in section IV as "FFLO-like", or 
quasi-one-dimensional states show a behavior 



a 
L 



X 



(82) 
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which is in agreement with this possibihty. The numer- 
ical value of the constant x is close to 2-^/2, which cor- 
responds to a = 27r^o and to the lattice constant n/q of 
the LO state. Thus, the LO state may be identified als 
the limiting case of the quasi-one-dimensional states of 
section IV for large n; the distance of the FFLO-lines is 
essentially independent of n, while the periodicity length 
6 sin a in the direction perpendicular to the lines tends 
to infinity (like hsma = rnr/q) for n ^ oo. The one- 
dimensional FFLO unit cell is a substructure that de- 
velops inside the diverging unit cell of the paramagnetic 
vortex states. 

To complete the description of the transition to the 
LO state, the above lattice structure may be used in 
Eq. (16) to perform the limit n ^ oo of the order param- 
eter expansion A„. We consider a 2D sample of finite 
area Fp, which contains NaNi, "small" unit cells of size 
Fc = ah sin a. The total area is given by Fp = LaL^ sin a 
with La = Natt and Lt = Ni,b. For different n the size 
and shape of Fc may change while Fp remains, of course, 
unchanged. Adopting the above model for the behavior 
of the unit cell as a function of n, we have n-independent 
numbers Na and a, while b — 6*^"^ increases linearly with 

(n) 

n and A^^ = N^^ decreases consequently according to 
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(83) 



Thus, a largest possible Landau number n = rig exists, 
which corresponds to 6^""=^ = Lf, (or N^"'"^ = 1) and is 
given by 



2Li) sin a 



(84) 



This cutoff Tic agrees exactly, in the present model, with 
the number of n = 1 unit cells fitting into a length Lb- 
As an additional consequence of the finite area of the 
sample, only a finite number of terms occur in the sum 
over m in Eq. (16). This number is fixed by the condition 
that the "center positions" y„i = m6sina = rmrL'^/a, lie 
inside the sample^^. This leads to the condition 



aLb sin a aLb sin a 

<m< 
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(85) 



which is in the limit n = ric only fulfilled for m = 0. 
Using the asymptotic expansion^^ of the Laguerre poly- 
nomial Ln, for large and even n = 2j, and taking into 
account only the term with m = in the sum of Eq. (16), 
the order parameter takes the form 



AC2jDj cos 



(86) 



The amplitude A [sec Eq. (15)] is, in the present system 
of (ordinary) units, given by 



1/4 



The coefficient C2j 
by C2j 

takes the form 



is, in the limit n 



simply given 



(Fp)i/2 jggg Eq. (17)], and the coefficient Dj 



' (-l)VM2j-l)!' 

While these factors. A, Dj, C2j diverge for n ^ oo, if the 
sample dimensions approach infinity, all singularities can- 
cel if n is replaced by the cutoff ric, and one obtains the 
expected result, = cos qy, for the one-dimensional 

periodic order parameter structure in the purely param- 
agnetic limit. 

The transition to the two-dimensional (square and 
hexagonal) periodic states found by Shimahara'^^ is more 
involved than the transition to the LO state. Let us re- 
strict to the square state, which is the simplest of all 2D 
states, and is also most stable from a thermodynamic 
point of view. 

For the square state, which is a linear combination of 
two LO states with orthogonal wave vectors, one would 
expect a divergent behavior of both unit cell basis vec- 



tors of the type a « n 



1/2 



,1/2 



Consequently, 



choosing a square unit cell in the (exact) order parameter 

expansion Eq. (18), one would expect to find a substruc- 
ture which becomes increasingly similar, with increasing 
n, to the structure of Shimaharas square state (line-like 
order parameter zeros, in the form of two sets of orthog- 
onal straight lines and circles). Numerical calculations, 
performed in the range n < 40 are, however, not in agree- 
ment with this expectation. 

On the other hand, the mathematical limit of the order 
parameter (18) yields in fact a 2D state with the period- 
icity of the FFLO wave vector and square symmetry, as 
shown in appendix E for a simplified model. The ex- 
planation for this apparent contradiction is provided by 
the result [relation (E8) of appendix E], that the quan- 
tum number n for a square state must obey the condition 
n = ttN'^, where A^ is an integer. This is a general result, 
which has been derived using essentially only the behav- 
ior a « n^/^ for large n. The latter is a consequence of 
the flux quantization condition and the shape of the unit 
cell. 

Of course, the relation n = ttN"^ cannot be fulfilled 
exactly for finite numbers n, N (for a sample of finite 
extension) since tt is an irrational number. The proper 
meaning of this relation is, that the sequence of states 
with quantum numbers n = int{'KN'^), N = 1,2,... rep- 
resents a sequence of approximations (of increased qual- 
ity) to the square state. Thus, the square state is the 
limit of a sequence defined on a very small subset of the 
set of integer numbers. 

This explains, why no systematic development of the 
square state with increasing n has been observed in the 
numerical calculations. The largest quantum number in 
the considered range (n < 40), which fulfills the above 
condition is n = 28 (corresponding to A' = 3). The order 
parameter modulus for n = 28 is shown in Fig. 10. It 
reveals, in fact, a certain similarity to the structure of 
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FIG. 10: (Reduced image quality due to arXiv restrictions) 
Contour plot of the square of the order parameter modulus 
for Landau quantum number n = 28 and a unit cell with 
parameters a = b, a = n/2. 



the square state (at least more similarity than any other 
state in the considered range). The arrangement of iso- 
lated order parameter zeros in Fig. 10 shows a tendency 
towards the formation of line-like zeros. Clearly, an ex- 
tremely high n and an extremely sharp definition of the 
tilt angle would be required to produce a really good ap- 
proximation to the square state. The final conclusion of 
the present analysis for the square state, that extreme 
requirements with regard to the definition of the tilt an- 
gle must be fulfilled in order to produce it, will probably 
hold for all other 2D states as well. 

The above analysis of the formation of the FFLO 
state(s) as limit(s) of the paramagnetic vortex states for 
n oo has been based on relation (79). In addition, 
relation (79) allows for an intuitive understanding of the 
unusual phenomenon of Cooper pairing at higher n, en- 
countered in the present configuration. The choice n = 
for the ordinary vortex state - in the absence of param- 
agnetic pair-breaking - corresponds to the lowest energy 
the system can achieve for p = 0. For sufficient large 
and decreasing H±, the Landau level spacing becomes 
smaller than the kinetic energy and the system has to 
perform a quantum jump from the n = to the n — 1 
pairing state, in order to fulfill the requirement of given 
energy as close as possible within the available range of 
discrete states [Inserting n = 1 in Eq. (79) determines 
the angle 6i as given by Eq. (1)]. For the same reason, a 
series of successive transitions to superconducting states 
of increasing n takes place with further decreasing H±, 
until the FFLO state is finally reached at B± = 0. The 
FFLO state for n ^ oo may obviously be considered as 
the continuum limit, or quasiclassical limit, of this se- 
ries of Cooper pair states, which starts with the ordinary 
vortex state at n = 0. 



VII. CONCLUSION 

The paramagnetic vortex states studied here, appear 
in a small interval of tilt angles close to the parallel ori- 
entation. A common feature of all of these states is a 
finite momentum of the superconducting pair wave func- 
tion, which is due to the large parallel component of the 
applied magnetic field. In these new superconducting 
states the Cooper pairs occupy quantized Landau levels 
with nonzero quantum numbers n. The number n in- 
creases with decreasing tilt angle and tends to infinity 
for the parallel orientation, where the FFLO state is re- 
alized. The unusual occupation of higher Landau levels 
may be understood in terms of the finite momentum of 
the Cooper pairs. 

The end points of the infinite series of Cooper-pair 
wave states occupying different n are the ordinary vor- 
tex state at n = and the FFLO state at n = oo. The 
dominant pair-breaking mechanism in the vortex state is 
the orbital effect, while Cooper pairs can only by bro- 
ken by means of the spin effect in the FFLO state. The 
equilibrium structure of the new states, which occupy the 
levels < 71 < cxD, is very different from the structure of 
the FFLO state(s), despite the fact, that the difference in 
tilt angles and phase boundaries may be small. Generally 
speaking, the equilibrium structures of the new states re- 
fiect the presence of both pair-breaking mechanisms; the 
fact that the local magnetic response may be diamag- 
netic or paramagnetic depending on the position in the 
unit cell may be understood in terms of this competition. 
A second unusual property, also closely related to the si- 
multaneous presence of both pair-breaking mechanisms, 
is the coexistence of vortices and antivortices in a single 
unit cell. 

The FFLO state has been predicted in 1964 and a large 
number of experimental and theoretical works dealing 
with this effect have been published since then. A def- 
inite experimental verification has not been achieved by 
now. However, recent experiments in the organic super- 
conductor K-{BEDT-TTF)2Cu{NCS)2 and other lay- 
ered materials'"'^'*^''*^''*^ revealed remarkable agreement'*'* 
with theory, both with regard to the angular- and the 
temperature-dependence of the upper critical field. In 
these phase boundary experiments, identification of the 
FFLO precursor states, studied in the present paper, 
seems possible if the tilt angle is defined with high pre- 
cision. To obtain a more direct evidence for all of these 
unconventional states, including the FFLO limit, other 
experiments, such as measurements of the local density 
of states by means of a scanning tunnelling microscope 
would be useful. 
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APPENDIX A: 



SYSTEM OF UNITS AND 
NOTATION 



In this appendix wc use primes to distinguish Eilcn- 
bergers dimensionless quantities, which wiU be used in 
sections III-V, from ordinary ones. The primes will be 
omitted in sections III-V. 

temperature: t = T/T^ 

length: r = r/Ro, Ro = TiVF/^TrkBTc = 0.882 ^o, 6 is 
the BCS coherence length. 
Fermi velocity: vp = vp/vp 
wave number: k = kRo 

Matsubara frequencies: = uJi/nkBTc = {21 + l)t 
order parameter: A = A/tt/cbTc 
magnetic field: H — H/Hq, where Hq = hc/2eRQ 
vector potential: A ~ A/Aq, where Aq — fi.c/2ei?o 
magnetic moment: /i = /i//io = T^kBTc/mVp, where 
fj,Q = irkBTc/ Hq. Note that the dimensionless magnetic 
moment agrees with the quasiclassical parameter. 
Gibbs free energy: G' = G/li-KkBTcf NpRl] 

Eilenbergers parameter R is related to the GL-parameter 
Kq of a clean superconductor according to the relation 

k= (j^C(3))'^''«o = 0.6837ko. 

The symbol k denotes a dimensionless, 2D unit vector. 
The Fermi-surface average of a fc— dependent quantity 
a{k) is denoted by a. For our cylindrical Fermi surface 
this average is simply an integral from to 27r over the 
azimuth angle ip. 

a = -J- (L 6^ka{k) = ^ [ dip a{k{ip)) . 
47r J 2-K Jq 

Finally, the symbol (a), defined by 



Ju 



unit cell 



denotes a spatial average of a quantity a(r) over a unit 
cell of area Fc- 



APPENDIX B: GAP CORRELATION 
FUNCTION 

It is convenient to express the gap correlation function, 
defined by Eq. (35) in terms of center of mass coordi- 
nates R = {r\ + r2) /2, r = r\ — r2, using the notation 
V'^^{R,r) = F(fi,r2). The function V^^{R,f) is in- 
variant under center of mass translations R => R+la+jb 
and may consequently be expanded in a Fourier series, 
using reciprocal lattice vectors Qij = IQi + jQ2, Uj = 
0, ±1, ±2, . . ., with basis vectors 



a 



1 




1 

sin a 



The Fourier coefficients of (i?, r) are denoted by 
V;j(r). The Fourier transform of Vij{r) with respect 

to r is denoted by v/^^p). 

Using the behavior of the gap A(r) under lattice trans- 
lations r f 4- rij, where = la + jb, the important 
relation 



y_;,,(f) = e'-'(^+^™^")yo,o(r + rj,i) 



(Bl) 



may be proven. This relation, first reported by Delrieu^^, 
shows that all Fourier coefficients are known if Vo,o is 
known. A similar relation holds for the Fourier transform 

Vffif) = e'P"^-.-'-''''(^+^™""Vo^J(p). 

The functions Vi_j and Vi^^\ which are most useful for 
the evaluation of the free energy, may be calculated by 
proceeding along the chain 

y^^(^,rl yo,o(r) V^jip) Vifip) Vtj{r), 

where an arrow denotes either calculation of a Fourier 
coefficient, or of a Fourier transform, or application of 
Delrieu's relation. 

Using the order parameter expansion (16) and per- 
forming the necessary manipulations, the result for Vj^J^ 
is given by 



/ 2 



(B2) 



where Qij,x, Qi,j,v are the x and y components respec- 
tively of the reciprocal lattice vector Qij. The final result 
for Vij is given by 



^,.(r) = (-l)'^dA„p). 



Ba 



01 



XL 



sma 



21 j cos a 



(B3) 



(B4) 



The usefulness of the gap correlation function for pair- 
wave states with with arbitrary n is essentially based on 
the translational invariance of the observable quantities 
\ib\ and B. 



APPENDIX C: THE GINZBURG-LANDAU 
LIMIT 

Let us first consider the upper critical field -ff^^ , which 
is determined by G^^^ = 0, for /z = 0, n = (if|| = 0) 
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and i — !■ 1. Solving this equation in this hmit, one finds, 
using ordinary units, 



i/f/ = 1.222— 



(CI) 



Eq. (CI) differs from the usual GL result by a factor of 
3/2. This discrepancy is due to our use of a cylindrical 
Fermi surface, instead of a spherical one, and can be elim- 
inated by replacing the GL parameter n by ?>k.qlI'2, (the 
quantities used in Eilenberger units are derived assuming 
a spherical Fermi surface). 

The magnetization relation (66) takes the following 
form for = 0, /i = 0, < ^ 1: 



B- H = AttM = 



H -H, 



c2 



2^2 /A J 



1 



(C2) 



The coefficient Aj_ in (C2)is given by (68). The fourth 
order free energy contribution (54) takes the form 



4 



4k2 



, (C3) 



where S'-^'^ = 7C(3)/8. The first sum in Eq. (C3) turns 
out to agree with Abrikosov's geometrical factor Pa, 



l,7n 



fl{xi,rn) = Pa, 



(C4) 



as discussed in more detail in KRS^. Performing again 
the above replacement of n one arrives at Abrikosovs 
well-known result 



dM 



(24^-l)/3^- 



(C5) 



Eq. (C4) remains also valid for n > 0. For the non- 
magnetic terms in Eq. (54), the Matsubara sum S^^^-^ 
may be considered as a low temperature correction to 
the GL term (C4). The GL- limit of the local magnetic 
field Bi± [see Eq. (63)] has also been calculated and has 
been found to obey the correct GL relation'^^ between 
magnetic field and square of order parameter. Here, the 
low-tcmperature corrections are contained in the Mat- 
subara sum sl 



(2) 



APPENDIX D: THE LIMIT OF THE ORDINARY 
VORTEX LATTICE 

It is of interest to investigate the limit of Eq. (54) cor- 
responding to the ordinary vortex lattice. We consider 
a situation without paramagnetic pair-breaking, i.e. set 
/Lt = 0, 9 = 7r/2, and ask for the equilibrium structure 
of the vortex lattice and the critical value of k separat- 
ing type I from type II superconductivity. To compare 
with the usual notation, we use here the same scaling 



K 2k/3 of the GL parameter as in appendix C. Fig- 
ure 1 1 shows the free energy G(4) 

as a function of a/L, a 
for K — 1.46 (k — 1.5) at t = 0.5. The flat minimum 
of G^''-' at a/L ~ 1.905, a = 60 indicates that the stable 
configuration is, as expected, a triangular vortex lattice. 
No other local minimum of the free energy exists. With 
decreasing k this minimum changes quickly into a maxi- 
mum; below K = 1.36 the free energy has no minimum at 
all which means that no spatially varying superconduct- 
ing state exists. The critical value of k = 1.36 separating 
type I from type II behavior at t = 0.5 agrees fairly well 
with the result of k = 1.25 obtained by Kramer^^ for 
the phase boundary between typ II and type II/ 1 behav- 
ior. For lower temperature the agreement is worse; at 
t = 0.2 the present theory gives k = 2.5 while Kramer's 
theory^^ gives k = 1.7. Recall that the error induced 
by the asymptotic approximation of subsection III D in- 
creases with decreasing temperature. 




FIG. 11: (Reduced image quality due to arXiv restrictions) 
Contour plot of the free energy G'-'*^ as a function of a/L 
and a without paramagnetic pair-breaking. Parameters are 
k — 1.5, t = 0.5, fj, = 0, G = 7r/2 . The minimal value 
G'"*' = 0.5067 is at a/L = 1.905, a = 60. 



APPENDIX E: THE SQUARE LIMIT FOR A 
MODEL ORDER PARAMETER 



The square of the the order parameter modulus, 
Eq. (18), for a square lattice may be written in the form 

|V„P(x,y,a) = Xi7,,,, (El) 

1,3 

Hi^, = (-l)'^e-f (''+j")L„(7r(/2 + j2)y^(i-+3v) _ 

We are interested in the limiting behavior of Eq. (El) for 
n — > oo, a « 77.2 — > oo_ In this limit, the quantity 27r/a 
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tends to zero and the double sum may be approximated 
by a double integral. An appropriate tool to perform 
such a calculation for infinite sums in a systematic way 
is Poissons summation formula. Using a two-dimensional 
version, which is derived in exactly the same way as for 
single sums, Eq. (El) may be written in the form 



(E2) 



\i^n\\x,y,a)={^yj2Y.f'^K UK 

rrix my 



where h{kx,x,ky,y,a) is a function representing Hij. 
The problem here is the factor (—1)'^ [see Eq. (El)], 
which must be represented by an infinite series of step 
functions^^. 

Since we arc more interested in the qiiestion if a limit 
with the correct periodicity and symmetry exists, than 
in the detailed functional form of this limit, we repre- 
sent the factor (—1)'-' approximately by the real part of 
e^xpia^ kxky / An , i.e. we use the function 



h{kx,x,ky,y,a) = Hij 



(E3) 



to represent Hij. Using this model, the absolute value 
of the r.h.s. of Eq. (E2) will be denoted by ^(a;, y, a) 
instead of |V'nP(a;, y, a). It takes the form 



(E4) 



S{x,y,a)=[^y\Y,Y.[<^K fdk'y- 

TYlx my 

h{k'x, X — ruxa, k'y,y — niya, a) | , 
where h{kx, x, ky,y, a) is given by 

h{kx,x,ky,y,a) = cos (^^^=oky^ e-m^{''v+''v) ■ 

Ln(^^{kl + kl)y<''--+''yy\ 

In order to perform the integrations, the relation 



' ■ m=0 ^ ' 

may be used to rewrite the Laguerre polynomial in the 

integrand as a sum of products depending on kx and 
ky separately. Then, the integration over kx may be 
performed^^ and, after a simple shift of the integration 
variable, a second relation'*^ 

(-l)"(")H2m(\/5)H2„-2m(V») 



may be used to calculate the sum over m. Performing 
the integration over kx one obtains the final result 



S{x,y,a) 
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277, 



-X e 



Hr, 



+ (_1)V^*^ 



(E5)| 



where the abbreviations x = x — rUxa, y = y — rUya have 
been used. 

We are interested in the limiting value of Eq. (E5) 
for n ^ 00. The asymptotic behavior of the Hermite 
polynomials^^ for large n implies 



(2n+l)i/2^a;_ 
a 



(E6) 



(2n + Xfl'^^m^ 



Since a « n^/^, the factor in front of x in Eq. (E6) re- 
mains finite for n ^ 00 and defines the FFLO wave vector 
q, i.e. 



(2n + 1)^2 



27r 



q, for n — > 00. 



(E7) 



Eq. (E7) implies a restriction on the possible quantum 
numbers n of a square FFLO state. The fact that an 
integer number A'' of wave lengths A = 27r/g' must fit into 
a length a, implies the condition 



n = ttN^ 



(E8) 



in the limit a — > 00. Condition (E8) is of a general na- 
ture and not a specific feature of our model. For quan- 
tum numbers n obeying Eq. (E8), all rux, m^— dependent 
phase factors in Eq. (E6) become multiples of 27r and may 
be omitted. The limit of S{x, y, a) for a —> 00 obtained 
in this way is well-defined, i.e. independent of any cutoff, 
and is given by 



lim S{x,y,a) | cos(ga;) cos(qy)|. 



(E9) 



A rotation of 7r/4 transforms Eq. (E9) into the more fa- 
miliar form'^^ \cos{q'x') + cos{q'y')\. The expected cor- 
rect result for Itpn]"^ is the square of the r.h.s. of Eq. (E9). 
Thus, the result of our model calculation differs from the 
exact result. A limiting state of the correct periodicity 
and symmetry has, however, been obtained. 
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